surfacePlotFill<-function (data, parameters, type = c("contour", "image", "persp"), 
    what = c("density", "uncertainty"), transformation = c("none", 
        "log", "sqrt"), grid = 50, nlevels = 20, scale = FALSE, 
    xlim = NULL, ylim = NULL, identify = FALSE, verbose = FALSE, 
    swapAxes = FALSE, ...) 
{
    grid1 <- function(n, range = c(0, 1), edge = TRUE) {
        if (any(n < 0 | round(n) != n)) 
            stop("n must be nonpositive and integer")
        G <- rep(0, n)
        if (edge) {
            G <- seq(from = min(range), to = max(range), by = abs(diff(range))/(n - 
                1))
        }
        else {
            lj <- abs(diff(range))
            incr <- lj/(2 * n)
            G <- seq(from = min(range) + incr, to = max(range) - 
                incr, by = 2 * incr)
        }
        G
    }
    grid2 <- function(x, y) {
        lx <- length(x)
        ly <- length(y)
        xy <- matrix(0, nrow = lx * ly, ncol = 2)
        l <- 0
        for (j in 1:ly) {
            for (i in 1:lx) {
                l <- l + 1
                xy[l, ] <- c(x[i], y[j])
            }
        }
        xy
    }
    data <- as.matrix(data)
    if (dim(data)[2] != 2) 
        stop("data must be two dimensional")
    densNuncer <- function(modelName, data, parameters) {
        if (is.null(parameters$variance$cholsigma)) {
            parameters$variance$cholsigma <- parameters$variance$sigma
            G <- dim(parameters$variance$sigma)[3]
            for (k in 1:G) parameters$variance$cholsigma[, , 
                k] <- chol(parameters$variance$sigma[, , k])
        }
        cden <- cdensVVV(data = data, parameters = parameters)
        z <- sweep(cden, MARGIN = 2, FUN = "*", STATS = parameters$pro)
        den <- apply(z, 1, sum)
        z <- sweep(z, MARGIN = 1, FUN = "/", STATS = den)
        data.frame(density = den, uncertainty = 1 - apply(z, 
            1, max))
    }
    pro <- parameters$pro
    mu <- parameters$mean
    sigma <- parameters$variance$sigma
    haveParams <- !is.null(mu) && !is.null(sigma) && !is.null(pro) && 
        !any(is.na(mu)) && !any(is.na(sigma)) && !(any(is.na(pro)))
    if (haveParams) {
        G <- ncol(mu)
        dimpar <- dim(sigma)
        if (length(dimpar) != 3) {
            haveParams <- FALSE
            warning("covariance must be a 3D matrix")
        }
        if (G != dimpar[3]) {
            haveParams <- FALSE
            warning("means and variance parameters are incompatible")
        }
        mu <- array(mu, c(2, G))
        sigma <- array(sigma, c(2, 2, G))
    }
    if (!haveParams) 
        stop("need parameters to compute density")
    if (swapAxes) {
        if (haveParams) {
            parameters$pro <- pro[2:1]
            parameters$mean <- mu[2:1, ]
            parameters$variance$sigma <- sigma[2:1, 2:1, ]
        }
        data <- data[, 2:1]
    }
    if (is.null(xlim)) 
        xlim <- range(data[, 1])
    if (is.null(ylim)) 
        ylim <- range(data[, 2])
    if (scale) {
        par(pty = "s")
        d <- diff(xlim) - diff(ylim)
        if (d > 0) {
            ylim <- c(ylim[1] - d/2, ylim[2] + d/2)
        }
        else {
            xlim <- c(xlim[1] + d/2, xlim[2] - d/2)
        }
    }
    if (is.null(dnames <- dimnames(data)[[2]])) 
        xlab <- ylab <- ""
    else {
        xlab <- dnames[1]
        ylab <- dnames[2]
    }
    if (length(grid) == 1) 
        grid <- c(grid, grid)
    x <- grid1(n = grid[1], range = xlim, edge = TRUE)
    y <- grid1(n = grid[2], range = ylim, edge = TRUE)
    xy <- grid2(x, y)
    if (verbose) 
        cat("\n computing density and uncertainty over grid ...\n")
    Z <- densNuncer(modelName = "VVV", data = xy, parameters = parameters)
    lx <- length(x)
    ly <- length(y)
    CI <- type
    DU <- what
    TRANS <- transformation
    if (length(CI) > 1) 
        CI <- CI[1]
    if (length(DU) > 1) 
        DU <- DU[1]
    if (length(TRANS) > 1) 
        TRANS <- TRANS[1]
    switch(EXPR = DU, density = {
        zz <- matrix(Z$density, lx, ly)
        title2 <- "Density"
    }, uncertainty = {
        zz <- matrix(Z$uncertainty, lx, ly)
        title2 <- "Uncertainty"
    }, stop("what improperly specified"))
    switch(EXPR = TRANS, none = {
        title1 <- ""
    }, log = {
        zz <- logb(zz)
        title1 <- "log"
    }, sqrt = {
        zz <- sqrt(zz)
        title1 <- "sqrt"
    }, stop("transformation improperly specified"))
    switch(EXPR = CI, contour = {
        title3 <- "Contour"
        filled.contour(x = x, y = y, z = zz, nlevels = nlevels, xlab = xlab,  color = cm.colors, asp = 1, ylab = ylab, main = "", ...)
    }, image = {
        title3 <- "Image"
        image(x = x, y = y, z = zz, xlab = xlab, ylab = ylab, 
            main = "", ...)
    }, persp = {
        title3 <- "Perspective"
        persp(x = x, y = y, z = zz, xlab = xlab, ylab = ylab, 
            theta = 60, phi = 30, expand = 0.6, main = "", ...)
    }, stop("type improperly specified"))
    if (identify) {
        TITLE <- paste(c(title1, title2, title3, "Plot"), collapse = " ")
        title(TITLE)
    }
    invisible(list(x = x, y = y, z = zz))
}



surfacePlotG<-function (data, parameters, type = c("contour", "image", "persp"), 
    what = c("density", "uncertainty"), transformation = c("none", 
        "log", "sqrt"), grid = 50, nlevels = 20, scale = FALSE, 
    xlim = NULL, ylim = NULL, identify = FALSE, verbose = FALSE, 
    swapAxes = FALSE, ...) 
{
    grid1 <- function(n, range = c(0, 1), edge = TRUE) {
        if (any(n < 0 | round(n) != n)) 
            stop("n must be nonpositive and integer")
        G <- rep(0, n)
        if (edge) {
            G <- seq(from = min(range), to = max(range), by = abs(diff(range))/(n - 
                1))
        }
        else {
            lj <- abs(diff(range))
            incr <- lj/(2 * n)
            G <- seq(from = min(range) + incr, to = max(range) - 
                incr, by = 2 * incr)
        }
        G
    }
    grid2 <- function(x, y) {
        lx <- length(x)
        ly <- length(y)
        xy <- matrix(0, nrow = lx * ly, ncol = 2)
        l <- 0
        for (j in 1:ly) {
            for (i in 1:lx) {
                l <- l + 1
                xy[l, ] <- c(x[i], y[j])
            }
        }
        xy
    }
    data <- as.matrix(data)
    if (dim(data)[2] != 2) 
        stop("data must be two dimensional")
    densNuncer <- function(modelName, data, parameters) {
        if (is.null(parameters$variance$cholsigma)) {
            parameters$variance$cholsigma <- parameters$variance$sigma
            G <- dim(parameters$variance$sigma)[3]
            for (k in 1:G) parameters$variance$cholsigma[, , 
                k] <- chol(parameters$variance$sigma[, , k])
        }
        cden <- cdensVVV(data = data, parameters = parameters)
        z <- sweep(cden, MARGIN = 2, FUN = "*", STATS = parameters$pro)
        den <- apply(z, 1, sum)
        z <- sweep(z, MARGIN = 1, FUN = "/", STATS = den)
        data.frame(density = den, uncertainty = 1 - apply(z, 
            1, max))
    }
    pro <- parameters$pro
    mu <- parameters$mean
    sigma <- parameters$variance$sigma
    haveParams <- !is.null(mu) && !is.null(sigma) && !is.null(pro) && 
        !any(is.na(mu)) && !any(is.na(sigma)) && !(any(is.na(pro)))
    if (haveParams) {
        G <- ncol(mu)
        dimpar <- dim(sigma)
        if (length(dimpar) != 3) {
            haveParams <- FALSE
            warning("covariance must be a 3D matrix")
        }
        if (G != dimpar[3]) {
            haveParams <- FALSE
            warning("means and variance parameters are incompatible")
        }
        mu <- array(mu, c(2, G))
        sigma <- array(sigma, c(2, 2, G))
    }
    if (!haveParams) 
        stop("need parameters to compute density")
    if (swapAxes) {
        if (haveParams) {
            parameters$pro <- pro[2:1]
            parameters$mean <- mu[2:1, ]
            parameters$variance$sigma <- sigma[2:1, 2:1, ]
        }
        data <- data[, 2:1]
    }
    if (is.null(xlim)) 
        xlim <- range(data[, 1])
    if (is.null(ylim)) 
        ylim <- range(data[, 2])
    if (scale) {
        par(pty = "s")
        d <- diff(xlim) - diff(ylim)
        if (d > 0) {
            ylim <- c(ylim[1] - d/2, ylim[2] + d/2)
        }
        else {
            xlim <- c(xlim[1] + d/2, xlim[2] - d/2)
        }
    }
    if (is.null(dnames <- dimnames(data)[[2]])) 
        xlab <- ylab <- ""
    else {
        xlab <- dnames[1]
        ylab <- dnames[2]
    }
    if (length(grid) == 1) 
        grid <- c(grid, grid)
    x <- grid1(n = grid[1], range = xlim, edge = TRUE)
    y <- grid1(n = grid[2], range = ylim, edge = TRUE)
    xy <- grid2(x, y)
    if (verbose) 
        cat("\n computing density and uncertainty over grid ...\n")
    Z <- densNuncer(modelName = "VVV", data = xy, parameters = parameters)
    lx <- length(x)
    ly <- length(y)
    CI <- type
    DU <- what
    TRANS <- transformation
    if (length(CI) > 1) 
        CI <- CI[1]
    if (length(DU) > 1) 
        DU <- DU[1]
    if (length(TRANS) > 1) 
        TRANS <- TRANS[1]
    switch(EXPR = DU, density = {
        zz <- matrix(Z$density, lx, ly)
        title2 <- "Density"
    }, uncertainty = {
        zz <- matrix(Z$uncertainty, lx, ly)
        title2 <- "Uncertainty"
    }, stop("what improperly specified"))
    switch(EXPR = TRANS, none = {
        title1 <- ""
    }, log = {
        zz <- logb(zz)
        title1 <- "log"
    }, sqrt = {
        zz <- sqrt(zz)
        title1 <- "sqrt"
    }, stop("transformation improperly specified"))
    switch(EXPR = CI, contour = {
        title3 <- "Contour"
        contour(x = x, y = y, z = zz, nlevels = nlevels,, col=grey(.4), xlab = xlab, ylab = ylab, main = "", ...)
    }, image = {
        title3 <- "Image"
        image(x = x, y = y, z = zz, xlab = xlab, ylab = ylab, 
            main = "", ...)
    }, persp = {
        title3 <- "Perspective"
        persp(x = x, y = y, z = zz, xlab = xlab, ylab = ylab, 
            theta = 60, phi = 30, expand = 0.6, main = "", ...)
    }, stop("type improperly specified"))
    if (identify) {
        TITLE <- paste(c(title1, title2, title3, "Plot"), collapse = " ")
        title(TITLE)
    }
    invisible(list(x = x, y = y, z = zz))
}
